Figure 5
rstudio
September 10, 2024 21:47:35 CEST
Introduction
Load libraries
library(Seurat)
library(canceRbits)
library(dplyr)
library(patchwork)
library(DT)
library(SCpubr)
library(tibble)
library(dittoSeq)
library(scRepertoire)
library(reshape2)
library(viridis)
library(tidyr)
library(grDevices)
library(ggpubr)
library(ggplot2)
library(ggnewscale)
library(RColorBrewer)
library(scales)
library(enrichplot)
library(clusterProfiler)
library("org.Hs.eg.db")
library(DOSE)
library(msigdbr)
library(stringr)Load single cell RNA-Seq data
# Load seurat object containing single cell pre-processed and annotated data, in the metadata folder
srat <- readRDS(params$path_to_data)
meta <- srat@meta.data
meta$WHO <- "SD"
meta$WHO[meta$patient %in% c("NeoBCC007_post", "NeoBCC008_post", "NeoBCC012_post", "NeoBCC017_post")] <- "CR"
meta$WHO[meta$patient %in% c("NeoBCC004_post", "NeoBCC006_post", "NeoBCC010_post", "NeoBCC011_post")] <- "PR"
srat <- AddMetaData(srat, meta$WHO, col.name = "WHO")
srat$WHO <- factor(srat$WHO, levels = c("CR", "PR", "SD"))colors
colors_B <- c(
"34" = "#ae017e",
"7" = "#377eb8",
"21" = "#f781bf",
"38" = "#fed976",
"15" = "#a6cee3" ,
"22" = "#e31a1c" ,
"3" = "#9e9ac8",
"4" = "#fdb462" ,
"24" = "#b3de69"
)
colors_clonotype = c("Small (1e-04 < X <= 0.001)" = "#8c6bb1",
"Medium (0.001 < X <= 0.01)" = "#41b6c4",
"Large (0.01 < X <= 0.1)" = "#fec44f",
"Hyperexpanded (0.1 < X <= 1)" = "#ce1256"
)
colors_Ig = c("IGHA1" = "#addd8e",
"IGHA2" = "#238443",
"IGHM" = "#dd3497" ,
"IGHD" = "#41b6c4",
"IGHG1" = "#fec44f",
"IGHG2" = "#fe9929",
"IGHG3" = "#ec7014",
"IGHG4" = "#fee391"
)Figure 5a
s <- subset(srat, subset = anno_l1 %in% c("B cells","Plasma cells"))
s@meta.data$seurat_clusters <- factor(s@meta.data$seurat_clusters, levels=c("34","7", "21", "38",
"15", "22","3", "4", "24"))
p <- SCpubr::do_DimPlot(sample = s,
reduction = "umap",
label = TRUE,
repel = TRUE,
label.color = "black",
legend.position = "none",
group.by = "seurat_clusters",
font.size = 25,
pt.size = 0.5,
colors.use = colors_B) + theme_minimal() + NoLegend() + theme(text = element_text(size=20))
p
Figure 5b
Idents(s) <- s$seurat_clusters
genes <- list( "B" = c("MS4A1", "CD19"),
"Ag presentation" = c("HLA-DRA", "HLA-DRB1", "HLA-DPB1", "CD40"),
"Cytotox" = c("GZMK", "PRF1"),
"Mem" = "CD27",
"Plasma" = c("IGKC", "IGHG1", "CD38", "SDC1", "JCHAIN"),
"Cell state" = c("G2M.Score", "S.Score","percent.mito")
)
p <- SCpubr::do_DotPlot(sample = s,
features = genes,
font.size = 18,
legend.length = 12,
legend.type = "colorbar",
dot.scale = 8,
colors.use = c("#7fbc41","#b8e186", "#f7f7f7","#fde0ef", "#c51b7d"),
use_viridis = FALSE
)
p
Figure 5c
q1 <- SCpubr::do_BarPlot( s,
group.by = "seurat_clusters",
split.by = "WHO",
position = "fill",
flip = TRUE,
order=FALSE,
legend.position = "none",
font.size = 32 ,
colors.use = colors_B,
xlab = "",
ylab = ""
) + xlab("Response") + ylab("% of B and plasma cells")
q1
DT::datatable(q1$data,
caption = ("Figure 5c"),
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c( 'csv', 'excel')))q1 <- SCpubr::do_BarPlot( s,
group.by = "seurat_clusters",
split.by = "Pathological.Response",
position = "fill",
flip = TRUE,
order=FALSE,
legend.position = "none",
font.size = 32 ,
colors.use = colors_B,
xlab = "",
ylab = ""
) + xlab("Response") + ylab("% of B and plasma cells")
q1
Open and pre-process clonotyoe data
This is the code used to create the combined expression
object from the BCR-Seq data. Due to personal confidentiality, the
object cannot be provided to GEO, but raw data are available on request
on EGA.
As default, the chunk related to BCR-Seq are not run.
# Run from raw BCR data stored in path_to_data or load the combined BCR data in the following chunk
# loading the data
vdj <- list()
dirs <- list.dirs(path = path_to_data, full.names = TRUE, recursive = FALSE)
samples <- c("NeoBCC007_post",
"NeoBCC008_post",
"NeoBCC012_post",
"NeoBCC017_post",
"NeoBCC004_post",
"NeoBCC005_post",
"NeoBCC006_post",
"NeoBCC010_post",
"NeoBCC011_post",
"NeoBCC014_post",
"NeoBCC015_post",
"NeoBCC018_post")
for (i in (samples)){
dir_TCR <- paste0(path_to_data, i, "_BCR_VDJ")
vdj[[i]] <- read.csv(paste0(dir_TCR, "/filtered_contig_annotations.csv"))
}
head(vdj[[1]])
# combine the contigs
combined <- combineBCR(vdj,
samples = gsub("_pre|_post", "", samples),
removeNA = TRUE, removeMulti = TRUE)
merge_combined <- bind_rows(combined, .id = "sample")if(params$accessibility == "unlock"){
combined <- readRDS(file.path("../",params$path_to_BCR))
srat <- RenameCells(srat,new.names = paste0(srat$patient, "_",gsub("(_1|_2|_3|_4|_5|_6|_7|_8|_9)*", "",colnames(srat))))
srat <- RenameCells(srat,new.names = gsub("_post", "",colnames(srat)))
srat$sample <- gsub("_post", "",srat$patient)
seurat <- combineExpression(combined, srat,
cloneCall = "strict",
group.by = "sample", chain = "both",
proportion = TRUE,
filterNA = TRUE)
s <- subset(seurat, subset = anno_l1 %in% c("B cells","Plasma cells") )
s@meta.data$cloneType <- factor(s@meta.data$cloneType, levels = unique(s$cloneType))
} else{
print("Please request access to the BCR-Seq data.")
}Figure 5d
Done on the highly confident B and plasma cells with exactly 2 IGH and IGL chains
if(params$accessibility == "unlock"){
df <- s@meta.data
df <- df %>%
separate( CTgene, c("IGH", "IGL"), sep = "_") %>%
separate( IGH, c("V", "D", "J", "C"), "\\.")
s <- AddMetaData(s, df$C, col.name = "Ig_level")
sss <- subset(s, subset = Ig_level != "NA")
dim(sss)
sss@meta.data$cloneType <- factor(sss@meta.data$cloneType, levels = c("Hyperexpanded (0.1 < X <= 1)", "Large (0.01 < X <= 0.1)", "Medium (0.001 < X <= 0.01)", "Small (1e-04 < X <= 0.001)"))
# define cell group membership
Idents(sss) <- sss$cloneType
sss@meta.data$Ig_level <- factor(sss@meta.data$Ig_level, levels = c("IGHA1",
"IGHA2",
"IGHM",
"IGHD",
"IGHG1",
"IGHG2",
"IGHG3",
"IGHG4"))
q1 <- SCpubr::do_BarPlot( sss,
group.by = "Ig_level",
split.by = "anno_l1",
position = "fill",
flip = TRUE,
order=FALSE,
legend.position = "right",
font.size = 30 ,
colors.use = colors_Ig) +
xlab("") +
ylab("Frequency of Ig isotypes")
q1
DT::datatable(q1$data,
caption = ("Figure 5d"),
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c( 'csv', 'excel')))
q2 <- SCpubr::do_BarPlot( sss,
group.by = "Ig_level",
split.by = "cloneType",
position = "fill",
flip = TRUE,
order=FALSE,
legend.position = "right",
font.size = 30 ,
colors.use = colors_Ig ) +
xlab("") +
ylab("Frequency of Ig isotypes")
q2
p <- q1|q2
DT::datatable(q2$data,
caption = ("Figure 5d"),
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c( 'csv', 'excel')))
} else{
p <- "Please request access to the BCR-Seq data."
}
Figure 5f
if(params$accessibility == "unlock"){
Idents(s) <- s$seurat_clusters
s@meta.data$seurat_clusters <- factor(s@meta.data$seurat_clusters, levels=c("34","7", "21", "38",
"15", "22","3", "4", "24"))
s@meta.data$cloneType <- factor(s@meta.data$cloneType, levels = c("Hyperexpanded (0.1 < X <= 1)", "Large (0.01 < X <= 0.1)", "Medium (0.001 < X <= 0.01)", "Small (1e-04 < X <= 0.001)"))
# define cell group membership
Idents(s) <- s$cloneType
de_markers <- DElegate::FindAllMarkers2(s, replicate_column = "patient", method = "edger", min_fc = 0, min_rate = 0.1)
signature_h <- de_markers$feature[de_markers$group1=="Hyperexpanded (0.1 < X <= 1)" ]
background <- rownames(s)
# enricher for cnet plot
GO_H_gene_sets = msigdbr(species = "human", category = "H")
msigdbr_t2g = GO_H_gene_sets %>% dplyr::distinct(gs_name, gene_symbol) %>% as.data.frame()
ego_h <- enricher(gene = signature_h, universe = background, TERM2GENE = msigdbr_t2g)
b <- barplot(ego_h, title = "Hallmarks´enrichment of hyperexpanded clones")
tmp <- b$data
tmp$ID <- str_trunc(tmp$ID, 100, "right")
tmp$ID <- factor(tmp$ID, levels = tmp$ID[order(tmp$Count)])
p <- ggplot(tmp, aes(Count, ID)) +
geom_bar(stat = "identity", color="black", fill = colors_clonotype["Hyperexpanded (0.1 < X <= 1)"]) +
theme_classic()+
geom_text(
aes(label = (paste( "padj=", format(p.adjust,scientific = TRUE, digits = 2)))),
color = "black",
size = 6,
hjust=1,
position = position_dodge(0.5) ) +
theme(text = element_text(size=18)) +
ggtitle("Hyperexpanded clones")
p
DT::datatable(p$data,
caption = ("Figure 3I"),
extensions = 'Buttons',
options = list( dom = 'Bfrtip',
buttons = c( 'csv', 'excel')))
} else{
p <- "Please request access to the BCR-Seq data."
}
Session Info
## R version 4.3.0 (2023-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.4 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=en_US.UTF-8
## [6] LC_MESSAGES=en_US.UTF-8 LC_PAPER=en_US.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Vienna
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] stringr_1.5.1 msigdbr_7.5.1 DOSE_3.26.2 org.Hs.eg.db_3.17.0 AnnotationDbi_1.62.2 IRanges_2.34.1
## [7] S4Vectors_0.38.2 Biobase_2.60.0 BiocGenerics_0.46.0 clusterProfiler_4.8.3 enrichplot_1.20.3 scales_1.3.0
## [13] RColorBrewer_1.1-3 ggnewscale_0.4.10 tidyr_1.3.1 scRepertoire_1.10.1 dittoSeq_1.12.2 canceRbits_0.1.6
## [19] ggpubr_0.6.0.999 ggplot2_3.5.1 viridis_0.6.5 viridisLite_0.4.2 reshape2_1.4.4 tibble_3.2.1
## [25] SCpubr_2.0.2 DT_0.32 patchwork_1.2.0 dplyr_1.1.4 Seurat_5.0.3 SeuratObject_5.0.1
## [31] sp_2.1-3
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.4 matrixStats_1.2.0 spatstat.sparse_3.0-3 bitops_1.0-7 HDO.db_0.99.1
## [6] httr_1.4.7 doParallel_1.0.17 tools_4.3.0 sctransform_0.4.1 backports_1.4.1
## [11] utf8_1.2.4 R6_2.5.1 vegan_2.6-4 lazyeval_0.2.2 uwot_0.1.16
## [16] mgcv_1.9-1 permute_0.9-7 withr_3.0.0 gridExtra_2.3 progressr_0.14.0
## [21] cli_3.6.2 spatstat.explore_3.2-7 fastDummies_1.7.3 scatterpie_0.2.1 labeling_0.4.3
## [26] sass_0.4.9 spatstat.data_3.0-4 ggridges_0.5.6 pbapply_1.7-2 yulab.utils_0.1.4
## [31] gson_0.1.0 stringdist_0.9.12 parallelly_1.37.1 limma_3.56.2 RSQLite_2.3.5
## [36] VGAM_1.1-10 rstudioapi_0.16.0 generics_0.1.3 gridGraphics_0.5-1 ica_1.0-3
## [41] spatstat.random_3.2-3 crosstalk_1.2.1 car_3.1-2 GO.db_3.17.0 Matrix_1.6-5
## [46] fansi_1.0.6 abind_1.4-5 lifecycle_1.0.4 edgeR_3.42.4 yaml_2.3.8
## [51] carData_3.0-5 SummarizedExperiment_1.30.2 qvalue_2.32.0 Rtsne_0.17 blob_1.2.4
## [56] grid_4.3.0 promises_1.2.1 crayon_1.5.2 miniUI_0.1.1.1 lattice_0.22-6
## [61] cowplot_1.1.3 KEGGREST_1.40.1 pillar_1.9.0 knitr_1.45 fgsea_1.26.0
## [66] GenomicRanges_1.52.1 future.apply_1.11.1 codetools_0.2-19 fastmatch_1.1-4 leiden_0.4.3.1
## [71] glue_1.7.0 downloader_0.4 ggfun_0.1.5 data.table_1.15.2 treeio_1.24.3
## [76] vctrs_0.6.5 png_0.1-8 spam_2.10-0 gtable_0.3.5 assertthat_0.2.1
## [81] cachem_1.1.0 xfun_0.43 S4Arrays_1.0.6 mime_0.12 tidygraph_1.3.1
## [86] survival_3.5-8 DElegate_1.2.1 SingleCellExperiment_1.22.0 pheatmap_1.0.12 iterators_1.0.14
## [91] fitdistrplus_1.1-11 ROCR_1.0-11 nlme_3.1-164 ggtree_3.13.0.001 bit64_4.0.5
## [96] RcppAnnoy_0.0.22 evd_2.3-6.1 GenomeInfoDb_1.36.4 bslib_0.6.2 irlba_2.3.5.1
## [101] KernSmooth_2.23-22 DBI_1.2.2 colorspace_2.1-0 tidyselect_1.2.1 bit_4.0.5
## [106] compiler_4.3.0 SparseM_1.81 DelayedArray_0.26.7 plotly_4.10.4 shadowtext_0.1.3
## [111] lmtest_0.9-40 digest_0.6.35 goftest_1.2-3 spatstat.utils_3.0-4 rmarkdown_2.26
## [116] XVector_0.40.0 htmltools_0.5.8 pkgconfig_2.0.3 sparseMatrixStats_1.12.2 MatrixGenerics_1.12.3
## [121] highr_0.10 fastmap_1.2.0 rlang_1.1.4 htmlwidgets_1.6.4 shiny_1.8.1
## [126] farver_2.1.2 jquerylib_0.1.4 zoo_1.8-12 jsonlite_1.8.8 BiocParallel_1.34.2
## [131] GOSemSim_2.26.1 RCurl_1.98-1.14 magrittr_2.0.3 GenomeInfoDbData_1.2.10 ggplotify_0.1.2
## [136] dotCall64_1.1-1 munsell_0.5.1 Rcpp_1.0.12 evmix_2.12 babelgene_22.9
## [141] ape_5.8 reticulate_1.35.0 truncdist_1.0-2 stringi_1.8.4 ggalluvial_0.12.5
## [146] ggraph_2.2.1 zlibbioc_1.46.0 MASS_7.3-60.0.1 plyr_1.8.9 parallel_4.3.0
## [151] listenv_0.9.1 ggrepel_0.9.5 forcats_1.0.0 deldir_2.0-4 Biostrings_2.68.1
## [156] graphlayouts_1.1.1 splines_4.3.0 tensor_1.5 locfit_1.5-9.9 igraph_2.0.3
## [161] spatstat.geom_3.2-9 cubature_2.1.0 ggsignif_0.6.4 RcppHNSW_0.6.0 evaluate_0.23
## [166] foreach_1.5.2 tweenr_2.0.3 httpuv_1.6.15 RANN_2.6.1 purrr_1.0.2
## [171] polyclip_1.10-6 future_1.33.2 scattermore_1.2 ggforce_0.4.2 broom_1.0.5
## [176] xtable_1.8-4 tidytree_0.4.6 RSpectra_0.16-1 rstatix_0.7.2 later_1.3.2
## [181] gsl_2.1-8 aplot_0.2.3 memoise_2.0.1 cluster_2.1.6 powerTCR_1.20.0
## [186] globals_0.16.3